Nonlinear evolution of surface morphology in InAs/AlAs superlattices via surface 

diffusion 
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Continuum simulations of self-organized lateral compositional modulation growth in InAs/AlAs 
short-period superlattices on InP substrate are presented. Results of the simulations correspond 
quantitatively to the results of synchrotron x-ray diffraction experiments. The time evolution of 
the compositional modulation during epitaxial growth can be explained only including a nonlinear 
dependence of the elastic energy of the growing epitaxial layer on its thickness. From the fit of the 
experimental data to the growth simulations we have determined the parameters of this nonlinear 
dependence, ft was found that the modulation amplitude don't depend on the values of the surface 
diffusion constants of particular elements. 



The process of self-organization during the growth of 
semiconductor epitaxial nanostructures is described us- 
ing two different models If there is a high density 
of monolayer steps on the vicinal surface (the crystallo- 
graphic miscut angle is larger than approx. 1°), a step- 
bunching instability occurs 2], but if the density of the 
the monolayer steps is low, a self-organized growth of 
two-dimensional or three-dimensional islands takes place. 
The latter process occurs, if the reduction of the strain 
energy due to an elastic relaxation of internal stresses 
in the islands outweighs the corresponding increase of 
the surface energy (morphological Asaro-Tiller-Grinfeld 
(ATG) instability 0,0, [|). 

The cited papers analyzed the self-organization process 
in a linearized approach from which a critical wavelength 
of the surface corrugation follows as function of material 
parameters. The exact nonlinear equation of the sur- 
face evolution was studied by Yang and Srolovitz Q and 
Spencer and Meiron [?[ for the case of a semi-infinite sub- 
strate. It was found that the shape of an initially har- 
monic surface waviness changes and a sequence of deep 
cusps is created. This behavior was observed using scan- 
ning electron microscopy (see, e.g., ja|). 

The physical properties of very thin layers (down to 
few monolayers) differ from the properties of the bulk. 
This difference leads to the creation of a stable two- 
dimensional layer at the surface (wetting layer) in the 
first stage of the Stranski-Krastanov growth mode. The 
occurrence of this so-called "wetting-effect" can be ex- 
plained by a nonlinear dependence of the elastic energy 
density on the layer thickness 0, [lj| • Simulations showed 
that the wetting-effect suppresses the growth of the cusps 
and subsequently it leads to the formation of surface is- 
lands [111]. These islands are unstable and coalesce (the 



Ostwald ripening 12 1) (l3| . However, numerical growth 
simulations indicate that an anisotropy of the surface en- 
ergy limits the ripening process and causes the creation 
of a nearly homogeneous array of islands (see Ref. 
among others, and the citations therein). 



The evolution of the surface morphology of multilay- 
ers has been studied only in a linearized approach so far 
[TBI From this approach, an unlimited growth of 

the modulation amplitude follows, which does not corre- 
spond to the experimentally observed stabilization of the 
modulation amplitude during the growth. The aim of 
this paper is to describe this stabilizing effect using the 
exact nonlinear equation of growth including the wetting- 
effect. We have simulated the time evolution of the spon- 
taneous lateral modulation of layer thicknesses in short- 
period semiconductor superlattices and we have found 
that the evolution of the modulation amplitude quan- 
titatively corresponds to the results obtained by x-ray 
scattering measurements. 

We have studied a series of InAs/AlAs superlattices 
grown by molecular beam epitaxy on InP (001) sub- 
strates. The substrate was covered by a 100 nm thick 
In 2 ,Al 1 _ ;r As buffer layer with the same chemical composi- 
tion as the average composition of the InAs/AlAs stack. 
The nominal thicknesses of both InAs and AlAs layers 
were 1.9 monolayers (mL) in all samples. The samples 
were prepared in a series with 2, 5, 10, and 20 super- 
lattice periods. The growth temperature was 530°C and 
the growth rate 0.5 mL/s. The details of the growth can 
be found elsewhere (l7| . 

X-ray grazing-incidence diffraction (GID) measure- 
ments were carried out at beamline ID01 at European 
Synchrotron Radiation Facility (ESRF) in Grenoble us- 
ing the x-ray wavelength 1.54 A. We have measured the 
diffusely scattered intensity distributions around 400 and 
040 reciprocal lattice points in the q x q y plane parallel 
to the sample surface. Two first-order lateral satellite 
maxima were observed at samples with 5 and more su- 
perlattice periods. An example of this experimental in- 
tensity map is plotted in Fig. [1] The distance of the 
lateral satellites from the specular crystal-truncation rod 
at q Xi y = is inversely proportional to the lateral modu- 
lation period; the period remains constant for all samples 
and was determined to be (267 ± 15) A. From the posi- 
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FIG. 1: Diffusely scattered intensity measured around the 
400 reciprocal lattice point on the sample with 20 superlattice 
periods. The contour step is 10°. 25. The diffraction vector is 
along the q x axis. The arrows show lateral intensity maxima 
mentioned in the text. 
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FIG. 2: The dependence of ei on the number of superlattice 
periods; this dependence describes the time evolution of the 
modulation during the growth. The circles with error bars 
are the experimental points obtained from the x-ray data, 
the full line represents the simulations. The dashed line is 
the evolution of ei calculated in the linearized case [TEL [Tsl] . 



tion of the satellite maxima in the q x % plane we deter- 
mined the modulation directions close to [310] and [130]. 
The amplitudes of the satellites increase with the growing 
number of the superlattice periods, while their widths de- 
crease. This indicates that during the epitaxial growth 
the periodicity of the lateral modulation improves and 
its amplitude increases. The details of the experimen- 
tal setup and results are described in the previous paper 
[HI], where we have determined the time- dependence of 
the modulation amplitude from the x-ray data. 

The diffuse scattered intensity can be calculated using 
the formula 



idiff(q) = A 



J d 3 rxh(r)e~*- u M 



(1) 



where A is a constant, \h is the crystal polarizability, 
h is the diffraction vector, and u is the displacement 
vector. Using the procedure described in the previous 
work [18|, we have extracted the correlation function 
e(x — x') = ((c(x) — c )(c(x') — c )}, from the measured 
data, where c(x) is the local InAs concentration averaged 
along the growth direction z and Co is the average InAs 
concentration. The dependence of the first coefficient £i 
of the Fourier series of e obtained from the experimental 
data is plotted in Fig. [2| this coefficient corresponds to 
the modulation amplitude. 

The evolution of the surface described by the function 
z = h(x, t) is driven by the surface diffusion and it can 
be described by the equation [lj| 



dh{x,t) D s 0fl 



Of 



V 2 ^i(x, t) + F + r/(x, i), 



(2) 



where D s the surface diffusion coefficient, 6 is number of 
atoms per unit area on the surface, Q, is atomic volume, 
kg is the Boltzmann constant, T is the temperature, jj, 



is the chemical potential on the surface, F is the depo- 
sition rate, and r\ is the deposition and diffusion random 
noise. The model is assumed to be independent on the 
third coordinate y and hence only one-dimensional sur- 
face modulation in a two-dimensional (x, z) space can be 
simulated. This approach is well justified for quantum 
wires or for strongly elongated quantum dots. For the 
quantum dots of a more symmetric shape this simulation 
can give limited information only. This effect is discussed 
below. 

The chemical potential fj, can be expressed as [2(| 



H(x,t) = -CjHmtjktlm 



-h(x) 



dh 



(3) 



where is the chemical potential of an ideally flat un- 
strained surface, 7 is the surface tension, k is the surface 
curvature, Cjkim are the components of the elastic tensor 
of the material, and is the strain tensor. The function 

fg? (h) describes a dependence of the elastic energy den- 
sity on the layer thickness h giving rise to a wetting-effect. 
For a Ge layer on Si, this function was approximated by 
the exponential function [20( 



f^\h) w 0.05 x E s (l - exp(-V/imi)), 



(4) 



where Eg is strain energy density in thick flat layer and 
h m \ is thickness of one monolayer. We have used an anal- 
ogous formula 



f^(h) = E w (l-eM-h/h w )), 



(5) 



where Eyy and hw are parameters depending on the lat- 
tice misfit and elastic constants of the layer. 

The strain energy was calculated by a direct solution of 
the linear elasticity equations using the boundary integral 
method. The method used is a multiple layer extension 
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FIG. 3: Simulated interface profiles in a In As/ Al As superlat- 
tice. Thick dashed lines denote the interfaces between InAs 
(under) and AlAs (above). Thin solid lines denote the inter- 
faces between AlAs and InAs. 



of the method in Ref . [|| for an isotropic continuum with 
periodic boundary conditions. The boundary conditions 
on the internal interfaces are described in Ref. 



2l|, the 



nearly singular integrals were calculated using method 
developed in Ref. [221 ]. 

The simulations have been performed with known ma- 
terial parameters. The surface energy 7, calculated from 
the first principles, was taken from the Ref. [23| as 
Urn -2 . The resulting structure of the interfaces inside 
the superlattice is shown in Fig. [31 From the simu- 
lations, the modulation period of 300 A follows, which 
is in a reasonable agreement with the observed value 
Lexp = (267 ± 15) A. It should be noted that the sim- 
ulated modulation period is affected by the size of the 
simulated region, since there can be only an integer num- 
ber of the waves in the simulated system of a given size. 
To eliminate the influence of the system size we have sim- 
ulated the growth of several systems of sizes 150, 225, 300 
and 400 nm. For various system sizes, the modulation pe- 
riods were always obtained in the interval (300 ± 20) A 
depending on the particular system size. The modulation 
amplitude is not affected by the system size. 

During the growth of first layers in the stack the modu- 
lation amplitude grows exponentially as predicted by the 
linearized theory [l5|, In the further growth stage 

however, the rate of the growth of the modulation am- 
plitude decreases (see Fig. [31 where the experimentally 
obtained values of the first Fourier coefficient s± of e are 
compared with the simulation results). 

The simulations show a good agreement with the 
experimental results in spite of the simplified one- 
dimensional model of the surface used. Transmission 



electron microscopy (TEM) on similar samples [17j re- 
vealed that the modulation is nearly one-dimensional in- 
deed, resulting in a quasiperiodic sequence of quantum 
wires; this fact explains why the one-dimensional model 
is sufficient for the simulation of the modulation kinet- 
ics. The TEM observations also demonstrated that if the 



average lattice parameter of the (relaxed) multilayer is 
larger than that of the InP buffer underneath (the actual 
multilayer structure is laterally compressed), the modu- 
lation direction is close to [100]; if the multilayer is lat- 
erally deformed in tension, the modulation direction is 
close to the crystallographic directions [310] and [130]. 
Of course, the one-dimensional model used here cannot 
predict the modulation direction. We ascribe the depen- 
dence of the modulation direction on the deformation 
sign to the anisotropic surface tension and anisotropy in 



elastic constants [24] . The degree of anisotropy of the 



surface tension is also affected by the actual strain in the 
layer [25| and this fact could therefore also explain dif- 
ferent modulation directions in the case of a tensile and 
compressive deformation of the multilayer. 

The continuum simulation also allows for the forma- 
tion of non-physical layers the thickness of which are 
fractional numbers of monolayers. However, our results 
based on a continuum approximation are qualitatively 
similar to the those obtained using an atomistic model 
and a monolayer step corrugation [2[ . 

The observed and predicted modulation periods 
roughl y co rrespond to the period given by the linearized 
theory [H HH, which prediction is L w 200 A. Accord- 
ing to the Ref. [2(| the surface diffusion of In is about 
50 times faster than Al, although an exact value of the 
surface diffusivity of In is not known. In [53] the sur- 
face diffusion constant of Al at 530° C was found to be 
1.5 x 10 cm 2 s _1 . The deposition flux of As atoms is 
higher than the flux of In and Al atoms at usual MBE 
conditions [28j, therefore only diffusivities of Al and In 
atoms play role. 

On the other hand, our simulations have shown that 
the values of the diffusion constants have nearly no in- 
fluence on the modulation amplitude, since the diffusion 
process is sufficiently quick and the growing surface is 
nearly in an equilibrium state. The diffusion rate how- 
ever affects the modulation period. In the case of very 
slow diffusion (of the order of 10 -10 cm 2 s _1 for Al), the 
growth of the larger ripples at the expense of smaller 
ones (the Ostwald ripening) does not take place and the 
modulation remains constant during the growth of the 
whole superlattice stack. If the diffusion is very fast (of 
the order of 10 _7 cm 2 s _1 for Al), the Ostwald ripening 
takes place during the growth of the first layer already, 
which leads to the creation of a smaller amount of larger, 
more distant dots, separated by larger flat areas of a thin 
wetting layer. The nucleation of the ripples on the sub- 
sequent interfaces is affected by the local distribution of 
lateral strains originated from the large buried ripples. 
Due to the elastic anisotropy, this distribution gives rise 
to local minima of the chemical potential at the rims of 
the buried ripples (two local minima for each ripple) so 
that the number of the ripples is duplicated. After the 
deposition of several periods, the ripples cover the whole 
interface again and the flat areas between the ripples dis- 
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appear. The resulting modulation period is approaching 
the period obtained for a slow diffusion again. 

In our simulations, we have achieved a good correspon- 
dence of both the modulation period and the time depen- 
dence of the modulation amplitude for any value of the 
diffusion constant of Al between 10 -10 and 10~ 7 cm 2 s -1 . 

The resulting interface morphology is substantially af- 
fected by the wetting-effect, i.e., by a non- linear depen- 
dence of the volume density of the elastic energy on the 
layer thickness. We have approximated this dependence 
by Eq. ((5J). The best correspondence of the measured 
and simulated modulation amplitudes was obtained for 
the values E w = 0.15 x E s and h w = 0.6 x h m \. We 
have also estimated these values by means of an atomistic 
simulation of the elastic energy density using the valence- 
field force method and the Keating model [29(. In these 
simulations we have neglected the surface relaxation and 
reconstruction and we have obtained the dependence of 
the density of the elastic energy on the thickness of a 
layer with a flat (001) surface. From the fit of this de- 
pendence with exponential formula in Eq. (|5|) we have 
obtained Ew = 0.10 x Eg and hw = 0.8 x h m \, which 
very well corresponds to the values above. 

The parameters Ew and hw affect the modulation am- 
plitude and they have no influence on the modulation 
period. In the first stage of the multilayer growth the 
modulation amplitude rapidly increases; this increase is 
slowed down after the growth of about 10 superlatticc 
periods. The parameter hw affects mainly the rate of 
the initial amplitude growth; this rate increase with de- 
creasing hw ■ The parameter Ew determines the slowing- 
down process: for larger values of Ew the slowdown of 
the amplitude growth is observed earlier than for smaller 
E W - 

In conclusion, we have simulated the multilayer growth 
using a non-linear continuum model. The simulation re- 
sults agree very well with experimental data obtained 
by x-ray scattering. From the simulations performed for 
various values of material parameters we have found that 
the wetting effect (the non-linear dependence of the elas- 
tic energy density on the layer thickness) has a crucial 
influence on the resulting interface morphology; from the 
fit of the experimental data with the simulations we have 
determined the parameters of this non-linear dependence 
and we have compared these values with atomistic simu- 
lations. 
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